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The random transitions of ion channels between conducting and non-conducting states generate 
a source of internal fluctuations in a neuron, known as channel noise. The standard method for 
modeling fluctuations in the states of ion channels uses continuous-time Markov chains nonlinearly 
coupled to a differential equation for voltage. Beginning with the work of Fox and Lu [1], there 
have been attempts to generate simpler models that use stochastic differential equation (SDEs) to 
approximate the stochastic spiking activity produced by Markov chain models. Recent numerical 
investigations, however, have raised doubts that SDE models can preserve the stochastic dynamics 
of Markov chain models [2-6] . 

We analyze three SDE models that have been proposed as approximations to the Markov chain 
model: one that describes the states of the ion channels and two that describe the states of the ion 
channel subunits. We show that the former channel-based approach can capture the distribution of 
channel noise and its effect on spiking in a Hodgkin-Huxley neuron model to a degree not previously 
demonstrated, but the latter two subunit-based approaches cannot. Our analysis provides intuitive 
and mathematical explanations for why this is the case: the temporal correlation in the channel 
noise is determined by the combinatorics of bundling subunits into channels, and the subunit- 
based approaches do not correctly account for this structure. Our study therefore confirms and 
elucidates the findings of previous numerical investigations of subunit-based SDE models. Moreover, 
it presents the first evidence that Markov chain models of the nonlinear, stochastic dynamics of 
neural membranes can be accurately approximated by SDEs, even with as few as 60 Na + and 18 
K + channels. This finding opens a door to future modeling work using SDE techniques to further 
illuminate the effects of ion channel fluctuations on electrically active cells. 



I. INTRODUCTION 



Hodgkin and Huxley's mathematical model of action 
potentials dynamics [7] is a cornerstone of computa- 
tional neuroscience. This system of equations provides a 
conductance-based framework for describing the dynam- 
ics of the membrane potential of a neuron. The essential 
features of the Hodgkin-Huxley (HH) model are quantita- 
tive descriptions of the permeability of a neuronal mem- 
brane to ion-specific currents (conductances) coupled to 
a current-balance equation that characterizes the voltage 
across a neural membrane. The physical basis for this 
empirical model is that conductances are determined by 
the proportion of ion channels in a conducting (i.e., an 
open state), and that the states of the ion channels are 
determined by the configuration of components of the ion 
channels, referred to as subunits, particles, or gates (sec 
[8, e.g.]). The original HH equations can be interpreted 
as a model of the mean behavior of the ion channels and 
their subunits. After the advent of techniques for record- 



ing the activity of individual ion channels [9] , it was dis- 
covered that ion channels can transition between open 
and closed states in an apparently random manner and 
that this can generate a significant source of noise to the 
ionic conductances [10]. 

This internal source of fluctuations is known as channel 
noise, and is distinguished from sources of noise that are 
external to the neuron such as synaptic noise [11, e.g.] 
or noise that is present in an applied stimulus. Chan- 
nel noise has important effects on neuronal dynamics and 
coding: it can alter the firing threshold [12-14], spike tim- 
ing [15, 16], interspike interval statistics [17], the amount 
of stochastic resonance [18, 19], and influence synaptic 
integration [20]. Channel noise can also contribute to 
the overall variability in the nervous system, which in 
turn may pose constraints on the fidelity of the motor 
and sensory systems of an animal [13, 21-24] and limit 
neuron miniaturization [25]. 

The classical HH formalism produces a deterministic 
description of neuronal dynamics, so alternative models 
have been proposed to account for channel noise. These 
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models assume that the activity of ion channels is gov- 
erned by random transitions among a number of possi- 
ble channel conformations, which leads to intrinsically 
stochastic models of neuronal dynamics. Although a va- 
riety of models of this type have been proposed, including 
those that capture fractal properties of patch-clamp data 
[26] and history-dependence in the activity of ion chan- 
nels [27], the most widely-used channel noise model is 
the Markov chain (MC) model. MC models assume that 
the state of an ion channel is described by a discrete- 
state, continuous-time Markov chain, where each state 
in the chain represents a particular configuration of the 
ion channel. The Markov property requires that a chan- 
nel's transition from one state to the next depends on its 
current state alone, thus the transition rates are deter- 
mined solely by the state of the channel and the voltage 
potential of the membrane. As a consequence, all chan- 
nels are coupled due to their common dependence on the 
membrane potential. For a recent review of MC models 
in computational neuroscience, see [28]. 

MC models are an invaluable tool for investigating the 
effects of channel noise on neural dynamics and coding, 
but these models are computationally expensive to sim- 
ulate and are difficult to analyze mathematically. As a 
result, there has been widespread interest in formulating 
stochastic differential equation (SDE) models of channel 
noise. This line of research was initiated by Fox and Lu 
[1, 29] and has been applied extensively in HH-type neu- 
ron models as well as models of calcium release from IP3 
receptors [30] (see [5] for a review of past applications of 
this approach). The SDE model that is most commonly 
used extends the original HH equations by including noise 
terms in the differential equations that describe the gat- 
ing variables. Computationally, this model can be orders 
of magnitude faster than MC models [2, 29], so it has 
often been used in place of, or as an approximation to, 
the MC model. Simulation studies have shown, however, 
that this SDE approach does not accurately replicate the 
stochastic response properties of the MC models [2-6] 
and it has been suggested that such models are inade- 
quate for simulations of channel noise [31] or must be 
modified to correctly reflect the stochastic properties of 
the MC models [5]. 

Despite recent indications that the commonly-used 
SDE model does not approximate the behavior of the 
MC model, there has been no definitive study detailing 
the cause of discrepancies between the MC and SDE ap- 
proaches. Moreover, other SDE models that have been 
proposed [1, 30] have never been tested to gauge whether 
there may be alternative, and more accurate, reduced 
models of channel noise. There are several possible rea- 
sons as to why an SDE model may not closely approx- 
imate a MC model. The system size expansion method 
for deriving an SDE model is an asymptotic result that 
is formally valid in the limit of a large number of chan- 
nels; it is possible that there are too few ion channels in 
a realistic model neuron for these approximate methods 
to be accurate. Another possible cause of the discrep- 



ancy between the two approaches is numerical error in 
the simulation algorithms [4]. Finally, it could be that 
the widely-used SDE models are formulated in a man- 
ner that neglects, or distorts, important dynamical and 
stochastic structure in the MC model. 

In this paper, we will demonstrate that the formulation 
of the SDE is critical for preserving the stochastic charac- 
ter of MC models. In Section II, we introduce three dif- 
ferent SDE models that have been proposed in the litera- 
ture. Among these, we distinguish between channel-based 
and subunit-based SDE models and provide an intuitive 
explanation for why the channel-based approach is the 
more appropriate SDE framework. We use a combination 
of mathematical analysis (Section III) and simulation re- 
sults (Section IV) to show that the MC model can, in 
principle, be well-approximated by a channel-based SDE 
model that was first introduced by Fox and Lu [1]. To 
our knowledge, ours is the first numerical implementa- 
tion of the channel-based Fox and Lu model [32] . Prior 
studies have provided numerical evidence that a widely- 
used subunit-based model does not accurately approx- 
imate the MC model [2-6], and our analysis confirms 
and elucidates these findings. We conclude that properly 
defining the structure and dynamics of ion channels is 
critical to formulating SDE models in a way that is con- 
sistent with MC models. We provide additional evidence 
for this conclusion by formulating reduced, quasistation- 
ary models based on our analytical results. Simulations 
of these models show that temporal correlation in the 
noise, which is shaped by the structure of ion channels, 
is critical for accurately approximating responses of the 
MC model. 



II. CONDUCTANCE MODELS BASED ON ION 
CHANNELS AND THEIR SUBUNITS 

We consider the HH model throughout this study, 
but our analysis is applicable to any conductance-based 
model with ion channels governed by linear, voltage- 
dependent kinetics. The membrane potential of an HH 
neuron is modeled as: 

C< ^t = ~ ffNa(F - ^ Na) - 9Kiy ~ Ek) ~ 9 ^ V ~ E ^ + I 

where C is the membrane capacitance, E^ ai E-^, and E^ 
are reversal potentials for Na + , K + , and leakage currents, 
respectively, and / is the applied current. Our central 
question is how to appropriately define the ion chan- 
nel conductances (<?Na for sodium and <?k for potassium) 
when one wants to include channel noise. Generally, one 
defines the K + conductance as <7k = </k/, where / is the 
fraction of K + channels that are open and <7k is the max- 
imal conductance per ion channel. Then the problem of 
appropriately reproducing K + channel behavior reduces 
to computing the evolution of /, the fraction of open 
channels. In the following, we will describe a number of 
methods for computing /. We outline the standard MC 
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model of ion channel kinetics [12, 33] and highlight how 
this approach relates to the classical (deterministic) HH 
model. We will then consider three distinct approaches 
for defining SDE models: two that were first proposed 
by Fox and Lu [1] and a variant suggested by Shuai and 
Jung [30]. 

Capturing the kinetics of a single subunit is the start- 
ing point for all of the models considered here. In the 
standard HH model, the K + channel has four indepen- 
dent identical subunits, traditionally given the symbol n, 
that must all be in an open state for the channel to be in 
the conducting state [1, 8]. The kinetics of a individual 
subunits are described by a two-state process: 



identical and independent subunits, its configuration can 
be modeled as a five-state Markov chain, where each state 
indicates the number of open subunits at a given instant 
in time: 



Aa n 3q„ 2a n a n 

fl 1 ^ 2 ^ 3 ^ 4. 

/3„ 2/3„ 3/9„ 4/3„ 



(4) 



The channel is said to be in the open or conducting state 
if all four subunits are open simultaneously. Let p be a 
column vector where the i th element represents the prob- 
ability at time t that a channel has i open subunits, then 
this probability distribution evolves in time according to 
the master equation 



Closed ^ Open, 

0n 



(2) 



where these subunits switch between open and closed 
states with the voltage-dependent transition rates [7]: 

0.01(10 - V) 



dp 

Itt 



(5) 



where the matrix A is: 



<*n(V) = 

P n (V) = 0.125 exp(-y/80). 

To simplify the notation, we will often omit the explicit 
dependence on V and write only a n and /?„. 

The Na + channels are modeled using two different sub- 
unit types, traditionally labeled m and h, each of which 
described by an open-closed kinetic scheme. The analysis 
of the two channel types is fundamentally the same, but 
entails significantly more notational complexity for the 
Na + channel. For conciseness, we will present detailed 
analysis of the K + channcl, but report results for both 
channel types. 

The analysis of the Na + channel can be carried out in 
an analogous fashion, but entails significantly more no- 
tational complexity, so we will present detailed analysis 
of the K + channel only. We report analytical and sim- 
ulation results for the both channels. In the remainder 
of this section, we review how this building block of a 
two-state subunit has been used to construct models of 
the K + conductance. 



A. Markov Chain Ion Channel Model 

The kinetic scheme in Eqn. 2 can be used to define 
a Markov chain that describes the behavior of a single 
subunit that randomly transitions between two states [8, 
e.g.]. If we let p sn ^ be the probability that the subunit is 
in the open state, then the evolution of this probability 
satisfies 
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= a n (l - p suh ) - /3 n p suh . 



(3) 



The conductance for a population of K + channels is 
determined by the proportion of the channels in the 
open state, <?k = <?k/, where <?k is the conductance 
per K + channel and / = N a /N is the fraction of open 
K + channels, with N a is the number of open channels and 
N is the total number of channels in the system. 



B. Deterministic Conductance Models 

If we consider an idealized neuron with an infinite num- 
ber of statistically-identical and independent channels, 
we can obtain a deterministic description of the fraction 
of open channels /. In this limit, the fraction of open 
channels is equivalent to the probability that any one 
channel will be open. In other words, Eqn. 5 also defines 
a deterministic model of conductance where <?k = ffKP4 
and p4 is given by the solution of the system of ordinary 
differential equations in Eqn. 5. 

At first glance, the deterministic definition of g K ap- 
pears to differ from that in the classical HH model, in 
which gK = g~YJ^ ■ As discussed in [8], however, these 
two models are equivalent: first, note that p su b in our 
notation can be identified with the gating variable n in 
the HH model because both satisfy the differential equa- 
tion (3) and both represent the proportion of subunits 
that are open. Next, observe that the entire system of 
differential equations in Eqn. 5 can be derived from the 
single HH gating variable by making the following sub- 
stitutions: 



This equation follows from Eqn. 2 and the fact that the 
probability that a subunit is closed is 1 — p S ub- Since 
the K + channel is assumed to consist of four statistically 



Pi 



: (!-«) 



4 -w, 



where i = 0, 1, 2, 3, or 4. 
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For instance, setting p4 = n 4 , we find: 
dp^ d 4 



dt dV 

= An 3 



, dn 
~dt 



= An 3 [a n (l - n) - /3 n n] 
= a n p 3 - Af3 n pi 

This equation is identical to the final row of Eqn. 5 and 
the remaining equations in that system (as well as those 
for Na + ) can be derived in a similar manner. 



C. SDE Conductance Models 

1. Channel SDE Model 

In the previous section we arrived at a determinis- 
tic model for K + conductance because we considered the 
case of infinitely many channels. If wc define the num- 



ber of K + channels to be finite, however, we can de- 
rive stochastic models using a system-size expansion [34] . 
This method was first applied to the HH model by Fox 
and Lu [1]. Following their notation, we define Xi to be 
the proportion of K + channels that have i open subunits. 
Since we are dealing with a finite population, the propor- 
tion of channels in a particular state, x i} is no longer a 
measure of probability pi. Rather, the number of open 
subunits fluctuates from realization to realization, which 
inevitably leads to a stochastic description of the chan- 
nel. The system size expansion provides a formal method 
for deriving a SDE model based on the master equation 
(5). Fox and Lu showed that the SDE for the K + channel 
is 



rfx 
Itt 



Ax + S£, 



(6) 



where x is a vector of the Xi, A is the matrix in Eqn. 5, 
£ is a vector of five independent Gaussian white noise 
processes with zero mean and unit variance, and S is the 
matrix square root of the diffusion matrix D, 



4a„xo+/5„ii — (4«„x +^„a;i) 

-(4a n x o +0„Xi) 4a n x + (3a n +l3 n )x 1 +2/3 n x 2 
-(3a n x 1 +2p n x 2 ) 





-(3a n x 1 +20 n x 2 ) 
3a n x 1 +2(a n +p n )x 2 +3/3 n x 3 
-(2a n x 2 +3/3 n x 3 ) 




2 a: 






-(2a n x 2 +3f3 n x 3 ) 
l x 2 + (a n +30 n )x 3 +4f3 rl x i 
-(a n x 3 +4l3 n x 4 ) 



-(a n x 3 +4f3 n x i ) 
a n x 3 +40 n X4 



(7) 



where N is the number of channels. To our knowledge, 
neither Fox and Lu nor other researchers have imple- 
mented this channel-based SDE model [32]. We will 
mathematically analyze it under voltage clamp condi- 
tions and perform numerical simulations to show that it 
accurately replicates the stochastic properties of the MC 
model. We refer to this model as the channel SDE model 
because the variable x is defined based on the states of 
the ion channels. 



2. Subunit SDE Models 

Unfortunately, the channel SDE model above does not 
preserve the dynamical structure of the classical HH 
equations: closed states are distinguishable and must 
be modeled, expanding the dimensionality of the sys- 
tem significantly. In an attempt to avoid this increase 
in complexity, one can apply the system size expansion 
procedure to the subunits rather than to the states of the 
channels. This leads to stochastic models that resemble 
the classical HH model, but include noise in the equations 
governing the subunit variables m, n, and h. We refer 
to such approaches as subunit SDE models. The subunit 
approach leads to the following SDE for the proportion 



of open subunits n [1, 30]: 
dn 

— =a n (l-n)- p n n + a n (V)t{t), (8) 

where the stochastic term is a Gaussian white noise 
process with zero mean and unit variance that is scaled 
in a voltage-dependent manner by a n (V), where 

^)= a " (1 -; )+ ^ n . o) 

The model for the channel population is then built 
from the subunit populations. Since each K + channel is 
composed of four statistically identical and independent 
subunits, Shuai and Jung proposed a model in which the 
proportion of open K + channels is defined as the prod- 
uct of four independent realizations (denoted rij) of so- 
lutions to the SDE in Eqn. 8 [30]. This defines the 
K + conductancc to be #k = 0k^i^2^3^4- We refer to 
this as the independent subunit model, InS. Shuai and 
Jung did not implement this method. Instead, they fol- 
lowed a method introduced by Fox and Lu [1] in which 
only one realization of a solution to Eqn. 8 is computed 
(denoted n), and this realization is raised to the fourth 
power. This defines the K + conductance to be <?k = 3k" 4 > 
built out of four identical subunit populations. We refer 
to this as the identical subunit model, IdS. In the limit 
of an infinite number of K + channels, both of the subunit 
models converge to the deterministic HH model. 
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D. The Distinction Between Subunit and Channel 
Models 

The fundamental difference between the channel SDE 
model in Section 11 C 1 and the subunit SDE models in 
Section II C 2 is that in the former, one first groups sub- 
units together to construct a channel and then defines 
the dynamics of the proportion of channels in each state. 
In the latter, one defines the dynamics of subunits first, 
before grouping the subunits together to compute the 
conductance of the channel. We note that the deter- 
ministic model in Section II B derived from the master 
equation is a channel-based approach while the classi- 
cal HH model is a subunit-based approach. Nonetheless, 
as discussed above, the two models are equivalent. It is 
tempting, therefore, to conclude that both the channel 
and subunit SDEs will also produce identical stochastic 
models. As we will show in the remainder of this study, 
these two approaches generate distinct stochastic pro- 
cesses: channel-centric SDE models can approximate the 
channel noise and spiking statistics of the MC model, but 
the subunit-based SDEs cannot. 

To gain some intuition for how the subunit and channel 
SDE approaches differ, consider the following example 
of a neuron with iV channels, where each channel con- 
sists of two statistically identical and independent sub- 
units. This configuration is illustrated in panel (a) of 
Fig. 1. The analysis can be extended to the four subunit 
K + channel, but for illustrative purposes we consider the 
simpler case of two subunits. At a given instant in time, 
define the state of the i th subunit in each class by the bi- 
nary random variables and Zi 2 . These variables take 
the value of one with probability p su b and are zero other- 
wise. The probability that the i th channel is open is de- 
termined by the probability that both subunits are open, 
p 2 sub - A channel-based approach defines the conductance 
from the proportion of open channels, so we average over 
all channels to obtain the proportion of open channels: 



(a) 



chan 



l N 



Zi,2 



(10) 



The Zi^Zi^ are themselves identical and independent bi- 
nary random variables that take the value one with prob- 
ability p 2 uh , thus the mean and variance of f c h&n are: 

E[f chan] 

Var[/ chan ] 



Psub 
1 



N 



Psub (! - Psub) 



The subunit approaches do not begin by grouping sub- 
units into predefined channels. Instead, one first com- 
putes the fraction of open subunits by averaging over each 
class of subunits, as shown in panel (b) of Fig. 1. The 
proportions of open subunits in the two subunit classes 
are: 



1 N 

•/sub ,j = Z i 
i=l 



where j = 1 or 2. 
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FIG. 1. Illustration of conceptual differences between 
channel-based and subunit-based models. In this example, 
each channel consists of two subunits: "A" and "O" ■ (a) 
In the channel-based model, subunits are first grouped to- 
gether to form channels (vertical rectangles) and the ionic 
conductance is determined by the fraction of channels in the 
conducting state, (b) In the independent subunits (IdS) ap- 
proach, the subunits are divided into two classes (horizontal 
rectangles) and the fraction of open subunits is computed by 
averaging over all subunits in each class. The proportions 
of open subunits in each class are then used to approximate 
the fraction of channels in the conducting state, (c) Iden- 
tical subunit (InS) models assume both subunits classes are 
identical. 



If we assume that each subunit class is independent and 
define the proportion of open channels / cnan to be the 
product of the / su b,i and / su b,2, then we can write the 
proportion of open channels as: 



fch&n 



1 



N 



i,3 = l 



The proportion of open channels under the subunit ap- 
proach is therefore an average of N 2 binary random vari- 
ables rather than TV random variables, as in Eqn. 10 for 
the channel approach. The probability that the product 
z i,\ z j,2 is equal to one is _p 2 ub , so the expected value of 
E[/ C han] = P 2 ub; identical to E[/ chan ]. The variance in the 
two models, however is different. To see this, write the 
variance of / c han as the sum of covariances: 



Var[/ ( 



chan I 



l 



N 



Y Cov(z M Z ij2 ,2: M 2: Zi 2) 
i,j,k,l=l 



This sum is over all of the TV 4 possible pairings of 
Zi,i, Zj,2, Zk,i, and z;^- To leading order in N, the domi- 
nant contribution to this sum is among pairings that have 
one index in common. There are N 2 possible pairs of 
and 3j 2: and for any given pair there are 2(N — 1) ways 
to choose the indices of Zk i and z; 2 such that i = k or 
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j = I. To leading order in N, therefore, there are 2N 3 of 
these terms, and the variance of / c han can thus be written 
as: 



Depending on whether 



is smaller or larger than 



Var[/ chan ] = 2pVar[/ chan ] 



(11) 



where p is the correlation coefficient between random 
variables of the form z^iZj^ and Zj^zip where either 
i = k or j = I. A straightforward calculation shows 
that p — i+p b b ■ Since p su b takes values in [0,1], p is 
bounded between zero and one half. Equation 11 implies, 
therefore, that the variance of the fraction of open chan- 
nels using the subunit-based averaging is always smaller 
than in the channel-based approach. Moreover, the vari- 
ance decreases with p su b- This implies that the subunit 
method underestimates the variance when p su b is small. 

An assumption in the above analysis is that there are 
two independent classes of subunits. This can be thought 
of as the analogue of the InS SDE model discussed above. 
Another variation of the subunit approach, used in the 
IdS SDE, is to assume that the two subunit classes are 
identical. This approach is illustrated in panel c of Fig. 
1. Since both subunit classes are identical, we drop the 
second subscript and define a random variable Zi that 
represents the state of the i th subunit in both classes. 
The proportion of open channels can then be rewritten 
as: 



/chan 



1 N 

= 



zl 



i=l 



1+Psub 

one fourth, therefore, the variance in the subunit model 
with identical subunit classes can be either smaller or 
larger than the variance of the MC model. These differ- 
ences are a direct consequence of how each approach ag- 
gregates the channels' subunits. Importantly, we observe 
that the differences between these approaches will persist 
for any finite number of channels. In the limit of infinitely 
many channels, the variance goes to zero, so all of the 
modeling approaches discussed here become equivalent. 
It is a straightforward exercise to extend this analysis 
to the case of four subunits (i.e., the K + channel), and 
a similar discrepancy between the channel and subunit 
approaches holds in that case. 

The three methods for grouping subunits that we have 
considered represent the three different approaches to 
performing a system-size expansion that have been pro- 
posed by Fox and Lu [1, 29] and Shuai and Jung [30]. 
These combinatorial arguments provide an intuitive un- 
derstanding for why the three SDE approaches that we 
are studying will lead to channel noise models with dif- 
ferent statistical properties. We now confirm this by di- 
rectly analyzing the SDE and MC models. 



III. VOLTAGE CLAMP ANALYSIS OF 
STOCHASTIC MODELS 

A. Stationary Distribution 



The expected value of /chan, which is given by the prob- 
ability that z\ — 1, is identical to E[/ cnan ]- The variance 
for this approach is: 



Var[/ ( 



chanj 



N 4 



N 

E 

;,j,k,i= 



Cav(ziZj,z k zi). 



As above, to leading order in N it suffices to consider the 
covariance of pairings that share one index in common. 
There is no distinction between the two subunit classes, 
so Cov(ziZj, ZiZk) = Cov(ziZj, z k Zi). There are therefore 
twice as many such terms as in the previous approach, 
and we find: 

Var[/ Chan ] = 4pVar[/ chan ] (12) 
with the same correlation coefficient p = , ^ 5 " b . 

r l+Psub 

This analysis illustrates how averaging across channels 
and averaging across subunits leads to fundamentally dif- 
ferent probabilistic descriptions of the proportion of open 
channels. In particular, since ,^ 3Ub < h, Eqn. 11 guar- 
antees that the variance of the proportion of open chan- 
nels given by the subunit model with two independent 
classes of subunits will never exceed the variance given 
by the MC model. Equation 12 shows that the variance 
in the subunit model with identical subunit classes will 
always be twice as large as the variance in the IdS model. 



We now analyze the stochastic character of each of 
these SDE models and compare them to the MC model. 
For the purposes of modeling voltage potential across the 
membrane, we are interested in the conductance of the 
ion channels so we seek to characterize the probability 
distribution of the fraction of open channels, which we 
denote /. To simplify the analysis, we will mimic the 
experimental technique of voltage clamp and perform our 
analysis while holding the membrane potential constant. 



1. Markov Chain Model 

In the MC model, each channel consists of four sub- 
units that transition between open and closed states. In 
voltage clamp this process is homogeneous in time so the 
stationary distribution for the number of open channels 
can be completely determined [34]. We are primarily 
interested in the stationary probability that all four sub- 
units are open because this is equivalent to the proba- 
bility that the channel itself is open. From Eqn. 3, the 
equilibrium value of p su b is - . The probability that 

the K + channel is open is therefore: 



Pchan — 



K + /3 n ) 4 ' 
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(a) Fraction of Open Na + Channels 

— O Markov chain 




-20 20 60 100 -20 20 60 100 

Membrane Potential (mV) 



(b) Fraction of Open K + Channels 




Membrane Potential (mV) 

FIG. 2. Analytical (lines) and numerical (symbols) results 
for means and standard deviations of the fractions of open 
channels in voltage clamp. The membrane area is 10 (im 2 
(600 Na + channels and 180 K + channels). The abscissa gives 
the voltage clamp value, (a) Results for Na + channels, (b) 
Results for K + channels. Analytical results for the channel 
SDE model are not shown because they are identical to those 
of the MC model. 

All K + channels are assumed to be statistically identical 
and independent, so the distribution of the total number 
of open K + channels at a given time is a binomial distri- 
bution with population parameter N (the total number 
of K + channels) and bias parameter p c han- To define the 
distribution of the fraction of open channels /, we rescale 
the binomial distribution by 1/N. The mean and the 
variance of this stationary distribution are: 

EmC [/] = Pchan — Mchan (13) 

tt r;i Pchan(l ~ Pchan) A 2 r-t »\ 

Var MC [/] = ^ = cr chan , (14) 

which we define as /i c han and crchan' respectively, as a 
shorthand. Note that these quantities are functions of 
V, even though we do not explicitly include this depen- 
dence in our notation. Equation 13 shows that the mean 
does not depend on the number of channels and Eqn. 14 
shows that the variance scales with 1/N. The mean and 
standard deviation of the fraction of open channels are 
shown by the solid black lines in Fig. 2. The mean and 
variance for the Na + channel can be computed in a simi- 
lar fashion and Fig. 2 includes those results. The results 
shown are for a membrane area of 10 /jm 2 , which corre- 
sponds to 180 K + and 600 Na + channels. 



2. Channel SDE Model 

To analyze the channel SDE model (Eqn. 6), we apply 
the simplification suggested by Fox and Lu: we set val- 
ues of the state variables in the diffusion matrix (Eqn. 7) 
to their mean equilibrium values. We refer to this ap- 
proximation as the equilibrium noise approximation and 
show in Appendix B that, for the voltage clamp case, 
it holds to 0(1/A^ 2 ). In principle, the values of each 
Xi should be confined to [0, 1] since they represent pro- 
portions of open channels, but to simplify the mathe- 
matical analysis we do not impose this condition. Un- 
der these simplifications, the SDE model is a multivari- 
ate Ornstein-Uhlenbeck (OU) process that by definition 
has a Gaussian-distributed stationary distribution. The 
mean and variance of the stationary process can be cal- 
culated directly using standard methods [34]. We find 
that the stationary distribution of the fraction of open 
K + channels is given by: 

1 (/->'ohan) 2 
P (f)= r-^- e -Ihan . (15) 
V 27rCT cha„ 

This is the expected result for a system size expansion 
since, in the limit of a large number of channels, the bino- 
mial distribution for the MC model can be approximated 
by a Gaussian distribution with mean /i c han and variance 
^chan- Any description of channel noise that aims to re- 
produce the behavior of the Markov state model should 
have this as its limiting distribution, so this result is our 
first confirmation that the channel SDE model provides 
an accurate approximation to the MC model. 



3. Subunit SDE Models 

To analyze the subunit-based SDE models, we apply 
the same approximations: we replace n in the equation 
for ex 2 (Eqn. 9) with its mean value at equilibrium. Thus, 
we replace with = N ^ n ^ s and do not restrict 

the values of n to the interval [0,1]. As was the case 
for the channel SDE, these approximations allow us to 
rewrite the SDE for n (Eqn. 8) as an OU process: 

dn 1 . 

— = — (/i sub - n) +a n £t, (16) 
at r n 

where r„ = an \p n and ^ sub = a "+p n ■ The stationary 
distribution of n (the fraction of open subunits) is there- 
fore Gaussian-distributed with mean fj, sn b and variance 

Vsub = N{*2+hn)* ■ Note that °s 2 ub is scaled b y 1/N; for 
simplicity we report analytical results to 0(1/N). 

In Section II C 2 we discussed two methods for defining 
the proportion of open K + channels based on the stochas- 
tic dynamics of the subunits. If we follow the approach 
of Shuai and Jung and combine four statistically identi- 
cal and independent solutions to Eqn. 16, the stationary 
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distribution for the proportion of open channels is de- 
fined by the product of four independent and identically 
distributed Gaussian random variables: 



pins(/) = J dm ...dm 



n 

3 = 1 



Unlike the channel SDE distribution in Eqn. 15, in the 
limit of a large number of channels, this distribution does 
not approach a Gaussian. This model, therefore, is fun- 
damentally incompatible with the MC model. Further- 
more, it is straightforward to compute the mean and 
variance of this distribution directly from the first two 
moments of the subunit distribution. We find: 



ElnS [/] = Msub 



Var InS [/] =4 M 6 -- 



Mchan 

sub^sub 



0(N~ 2 ) 



(17) 
(18) 



^ 2 



chan' 



We note that this leading order result for the moments 
can also be obtained following the combinatorial ap- 
proach outlined in Section II D. 

If we go further and assume that all four populations 
of subunits are identical and perfectly correlated (i.e., 
following the approximation proposed by Fox and Lu), 
then the stationary distribution for the fraction of open 
K + channels is given by the Gaussian distribution for a 
single subunit raised to the fourth power. The distribu- 
tion in this case has the closed form: 



Pids(f) = 



2 

sub 



/ 



-3/4, 



sub 



nt. 



This distribution also does not limit to a Gaussian and 
is fundamentally incompatible with the MC model. As 
above, we can compute the mean and variance of n . 



Ems [/] - Msub - 

Mohan 

Var IdS [/] = 16/4 



6M s Vsub + 0(N- 2 ) 
-OiN- 1 ), 

0(N~' 



2 

b^sub 



(19) 
(20) 



chan ' 



Equations 17 and 19, show that the mean fraction of 
open channels computed with the subunit approaches 
agrees with the MC model in the limit of a large number 
of channels. The variance, however, is poorly described. 
For instance, the ratio of the open channel variance of the 
MC model (Eqn. 14) to that of the IdS model (Eqn. 20) 
is: 

g chan = 1 + Msub + /4jb + flfub 

Var IdS (/) " 16 Ms 3 ub 

For subthreshold values of V, fi su b is small and therefore 
the IdS model drastically underestimates the magnitude 
of the channel noise. 



These analytical results for the two subunit SDE mod- 
els are plotted in Fig. 2 with dotted (independent sub- 
unit populations; InS) and dashed (identical subunit pop- 
ulations; IdS) lines. Note that the channel noise in the 
Na + channel is also underestimated by both subunit mod- 
els for V near the resting potential of zero millivolt. 



B. Autocorrelation in Voltage Clamp 

We now analyze temporal correlations in the propor- 
tion of open channels for a given voltage clamp level. 
As in prior sections, equations presented here depend on 
voltage potential V, but to simplify notation we do not 
explicitly indicate this. If we denote the time series of 
the proportion of open channels as f(t), then the auto- 
correlation function for f(t) is: 



R(t) 



E[/(t)/(0)]-E[/(0)]S 
Var[/(0)] 



We assume R(t) does not depend on the initial time since 
our analysis is restricted to the stationary distribution of 
open channels in voltage clamp. 



1. Markov Cham Model 

Let Ci(t) denote the state of the i th channel at time t, 
where Cj(i) = 1 indicates an open channel and Ci(t) = 
indicates that the channel is closed. The autocorrelation 
for the fraction of open channels then becomes 



R(t) = 



E 



^£?j=iCi(f)ci(0) 



II 2 

Mchan 



(21) 



chan 

E[ Ci (*) Ci (0)]-/4 an 



chan 



This simplification is possible because the MC model as- 
sumes that all channels are statistically identical and 
independent. The only unknown term in Eqn. 21 is 
E [cj(t)cj(0)]. Since c$ is a binary random variable, ex- 
pected value of Ci(t)ci(0) is equal to the probability that 
the channel is open at the initial time and is also open 
at the later time t. This probability can be determined 
by solving the master equation (5), which is possible in 
voltage clamp because this system of ordinary differential 
equations is a linear equation with constant coefficients. 
The probability that the channel is open is given by 
in Eqn. 5, so E [cj(i)cj(0)] is equal to the entry in the last 
row and the last column of the matrix exponential of the 
matrix in (5). We find: 



R(t) 



Aa? n e r n +6 a 2 i f3 n e ^ + 4q»/3^e ^ + fije ^ 

(2a n + (3 n )(2al + 2a n f3 n + ffi) ' 

(22) 

— T-5-- The solid black line in panel (b) 



where r„ = — nr 
of Fig. 3 shows this function for a voltage clamp value 
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FIG. 3. Analytical (lines) and numerical (symbols) results for 
autocorrelations R(t) of the fractions of open channels in volt- 
age clamp. The voltage clamp is set to mV and the mem- 
brane area is 10 /j,m 2 (600 Na + channels and 180 K + channels). 
(a) Results for Na + channels. (b) Results for K + channels. 
Analytical results for the channel SDE model are not shown 
because they are identical to those of the MC model. 



of V = mV. The same analysis was applied to the 
Na + channel, the result of which is also shown in panel 
(a) of Fig. 3. Note the difference in time scales on the 
ordinate in subpanels (a) and (b) . Temporal correlations 
in the fraction of open Na + channels decay rapidly within 
the first millisecond, whereas the temporal correlation in 
the fraction of open K + channels persists for nearly 10 ms. 

Equation 22 reveals an important feature of the MC 
models — the structure of the channel defines the tempo- 
ral profile of the channel noise statistics. In the case of 
the K + channel, the transitions between the five possible 
channel configurations induce correlations on four time 
scales, the first four multiples of l/r„. The Na + channel 
has eight possible channel configurations because it has 
three m subunits and one h subunit, so the autocorre- 
lation function for the fraction of open Na + channels has 
seven time scales. 



2. Channel SDE Model 

Using the same approximations that allow the channel 
models to be described as a multivariable OU process, 
we can compute the autocorrelation function for the pro- 
portion of open channels using standard methods [34]. 
We find that the autocorrelation function is identical to 
Eqn. 22 and thus do not include it in Fig. 3 



3. Subunit SDE Models 

We compute the autocorrelation for the InS and IdS 
models using the OU approximation of n in voltage 
clamp. We derive our results using the Ito calculus. The 
well-known solution to the OU process in Eqn. 16 at long 
times is [34, e.g]: 



l(t) = Msub + 



dW u . 



(23) 



To calculate the autocorrelation in the InS model, we 
take the expectation of the product of four independent 
solutions of the form of Eqn. 23 and normalize by the 
voltage clamp mean and variance shown in Eqns. 17 and 
18: 



■Rlns(*) 



E[ILf =1 m (t)«i(0)] - E InS (/) 2 



4 M ; 



sub £ 



Var InS (/) 
" + 6^ ub <7 



6 

sub 



6^; 



4 a 2 

sub sub 



+ <>(&). (24) 



For brevity, we omit terms of order higher than 1/N. 

For the IdS model, we instead take a single solution 
of the form of Eqn. 23 and raise it to the fourth power 
and normalize by the voltage clamp mean and variance 
shown in the Eqns. 19 and 20. The autocorrelation of 
this function is: 



E[n(t) 4 n(0) 4 ]-E IdS (/)- 

Var Ids (/) 
(2^ u b + 12^ ub a s 2 ub )e- T 



9/4 



ub^sub^ 



2/4b 



21fi suh a suh 



(25) 



Calculation of the higher order terms shows that the 
same exponential time scales (the first four multiples of 
l/r„) are present in R(t) for all of the models, but the 
coefficients arc different. In particular, Eqns. 24 and 25 
show that, to leading order in 1/N, the subunit models 
lack the faster time scales of the K + channel. As a result, 
temporal correlations in the subunit-based channel noise 
models persist longer than those of the MC model and 
the channel-based SDE model. Fig. 3, voltage clamped 
at V — mV, displays these differences in the autocorre- 
lation functions of the subunit models. The dotted line 
shows the result for the InS model (independent subunit 
populations) and the dashed line shows the result for 
the IdS model (identical subunit populations). The de- 
cay time in the autocorrelation of the proportion of open 
channels is longer for the subunit models for both the 
Na + channel shown in panel in (a) and the K + channel 
shown in panel (b). 



IV. NUMERICAL SIMULATIONS 

In this section we report results from numerical simu- 
lations of the MC model as well as the three SDE models 
analyzed above. We first verify the results of our analysis 
of voltage clamp statistics and then measure the statis- 
tics of interspike intervals in order to test how well the 
SDE models replicate the stochastic features of the MC 
model when voltage is allowed to evolve freely according 
to Eqn. 1. In all simulations we use the parameter val- 
ues listed in Table I. We perform simulations for three 
different membrane areas: 1, 10, and 100 /im 2 . The cor- 
responding channel counts are shown in Table II. 
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TABLE I. Parameter values from [7]. Note that resting po- 
tential has been shifted to mV. 



oy iiiuui 


T)p~h n 1 ti r»n 

i/L 11111 LlUll 


\/a nip iitii^q 


G 


Membrane capacitance 


1 fib J cm 


<?Na 


Maximal sodium conductance 


12U ma/cm 


m 


Maximal potassium conductance 


06 mS/cm 




Leak conductance 


0.3 mS/cm 2 




Sodium reversal potential 


115 mV 


Ek 


Potassium reversal potential 


-12 mV 


£l 


Leak reversal potential 


10.6 mV 


PNa 


Sodium channel density 


60 / fim 2 


PK 


Potassium channel density 


18 / fim 2 



TABLE II. Membrane areas and corresponding channel 
counts used in numerical simulations 



Membrane Area 


# of channels 


(/Mil 2 ) 


Na+ K+ 


1 


60 18 


10 


600 180 


100 


6000 1800 



A. Methods 

Sample Fortran code used to simulate the four stochas- 
tic models is available at the ModclDB website (acces- 
sion number 128502) [35]. All simulations used a time 
step of size 0.01 ms. We define a spike using two condi- 
tions: the membrane potential must exceed 60 mV and it 
must have remained below 60 mV for the previous 2 ms 
(approximately the width of a spike). To generate Gaus- 
sian random numbers, we first produced uniform random 
numbers using the Mersenne Twister algorithm [36] and 
then transformed these to Gaussian random numbers us- 
ing the Box-Muller method [37]. 

In Section IV C, we characterize the spiking response 
of the different models in response to stimuli of the form: 



I(t) — J D C + Inoise^ti 



(26) 



where £t is a Gaussian white noise process with zero 
mean and unit variance. This type of input is commonly 
used to characterize the response of stochastic Hodgkin- 
Huxley models [6, 17, e.g.]. The additive white noise 
term can be interpreted as a simplified method for repre- 
senting the combined effect of numerous synaptic inputs 
that neurons in cortex and other networks receive in vivo, 
see for instance [38] . We simulate spike trains for varying 
membrane area, DC input current Idc, and input noise 
Inoise- We report the mean and coefficient of variation 
(CV) for the first 2000 interspike intervals (ISI) for each 
spike train. 



1. Markov Chain Model 

The Markov chain describing each K + channel is shown 
in Eqn. 4. The Markov chain that governs the state 
of each Na + channel includes three m subunits and one 
h subunit and is therefore described by an eight state 
Markov chain: 



0,0 


3a m 


1,0 


2a m 
0m 


2,0 


a m 

Pm 


3,0 




0,1 


3a m 


cthltPh 
1,1 


2a m 


2,1 




3,1 


(27) 



The channel is in the conducting state when all three 
m subunits and the h subunit are open. The voltage- 
dependent transition rates for the m and h subunits are 
[7]: 



a m (V)=0.1 



25 -V 



exp( 



25-y 

10 



)-l 



l (V)=4exp(-g) 



a h {V)=0.m eX p(~) 



exp( 



3o-y 

10 



+ 1) 



The Markov chains in Eqns. 4 and 27 define the possi- 
ble states of each individual channel. Rather than simu- 
lating individual channels in the membrane patch, how- 
ever, it is more efficient to track the number of channels 
in each state using the Gillespie algorithm [39, 40]. At 
each time step, the fraction of open Na + channcls /n & and 
K + channels /k is computed and the voltage is updated 
using the forward Euler algorithm applied to Eqn. 1. 



2. Channel SDE Model (US) 

The channel SDE model is a system of twelve differen- 
tial equations derived by Fox and Lu [1]. In matrix form, 
it can be written as: 

dV 

C— = -m*V3i{V - E Na ) - g K x 4 (V - Ek) - g L (V - E L ) + 1 



rfx 
~dt 
dy 
dt 



A K x + 4a„x ei + S k £k 

^N a y + 3a m y oei + a h yoo^A + SW^Na- 



(28) 



The vector x is made up of entries Xi (i = 1,2,3,4) 
that represent the proportion of K + channels with i open 
n subunits. The entries of y are denoted [i = 
0,1,2,3 and j = 0,1) and represent the proportion of 
Na + channels with i open m subunits and j open h sub- 
units. The vectors ei are column vectors with a one in 
the i th entry and zero elsewhere. Following Fox and Lu, 
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we use the fact that X^=o Xi = 1 ana - Si=o Sj=o ^1? = 1 
to define xo and yoo- This allows us to reduce the dimen- 
sion of the system of SDEs from fourteen to twelve. We 
note that this reduction of dimension is exact, following 
from properties of the A and S matrices. We include 
the details of the matrices ^Na, ^4k, SW and 5k in the 
Appendix A. 



The values of the Xi and represent proportions of 
channels in a particular configurations so, in theory, they 
should lie within five-dimensional and eight-dimensional 
hypercubes bounded by the intervals [0,1]. Moreover, 
since the xi and each sum to one, the values of these 
variables should in fact lie on hyperplanes within these 
hypercubes. If, in the course of numerical simulations, 
there are excursions of x and y off of these high dimen- 
sional bounded surfaces, then the solution will lack bio- 
logical meaning because a value that represents the pro- 
portions of channels should not be negative or exceed 
one. A numerical difficulty also arises when these vari- 
ables do not lie on the proper bounded hyperplanes. To 
define Sk and S^a, one needs to take matrix square roots 
of diffusion matrices in every time step. If the x and y 
do not lie on the bounded surfaces, then the diffusion 
matrices, which depend on the values of x and y, will no 
longer be guaranteed to be positive semi-definite, which 
may make it impossible to compute real valued matrix 
square roots. 



In principle, it may be possible to incorporate a projec- 
tion or a reflection into the numerical method to ensure 
that x and y remain on these bounded, high-dimensional 
surfaces of admissible values. We demonstrate that a 
simpler approach in which the individual values of the Xi 
and are not confined within [0, 1], but rather are free 
to evolve without boundary conditions, gives an adequate 
numerical approximation to the intcrspike interval statis- 
tics of the MC model. With this simplification, there is 
no longer a guarantee that real matrix square roots of the 
diffusion matrices will exist, so we replace the values of 
Xi and ytj in the diffusion matrices with their equilibrium 
values. The validity of approximation is discussed in Ap- 
pendix B. After implementing the above simplifications, 
we solved the resulting system of SDEs using the Euler- 
Maruyama method [41]. Comparing our results with im- 
plementations that bound the SDE solutions would be 
an interesting subject for future work, but is beyond the 
aims of this paper. 



3. Subunit SDE Models 

The two subunit SDE models that we study are the 
independent subunit (InS) model: 

dV 

C— = -gN a m 1 m 2 m 3 h(V - E Na ) 



drrii 

~dT 
dh 

~dt 

drii 

~dt 



- 9Knin 2 n 3 n 4 (V - E K ) - g L (V - E L ) + I 

= a m (l - raj) - /3 m TOj + P m (V)€ mi (t), where i = 1,2,3 

= a h (l-h)-0 h h + v h (V)Z h (t) 

= a n (l - Hi) - f3 n rii + <r n {V)€ ni (t), where i = 1,2,3,4. 



and the identical subunit (IdS) model 

°1I = -9^™ 3 HV - #Na) - 9ku\V - E K ) - g L (V - E L ) + I 



dm 

~dt 
dh 

~dt 
dn 

~dl 



= a m (l - ra) - /3 m m + <r m (V)£ m (t) 
= a h (l-h)-p h h + a h (V)^ h (t) 
= a n (l - n) - (3 n n + <r n {V)£ n (t). 



The difference between these two models is that, in the 
former, we compute multiple independent realizations of 
the n and ra type subunits and the product of these 
terms enter into the equation for V whereas in the lat- 
ter, the approximation is made that all subunit classes 
are assumed to be perfectly correlated so only one SDE 
is solved for each subunit type and the solution is raised 
to the appropriate power (4 for n and 3 for ra) . The gat- 
ing variables represent proportions of open subunits so 
we enforce boundary conditions that prevent the values 
of the gating variables from exceeding one or becoming 
negative. 

We note that the form of a 2 x [x = in, h, or n) that 
we use is given in Eqn. 9. In particular, the noise terms 
depend on voltage as well as the subunit variables them- 
selves. We do not apply the equilibrium noise approxi- 
mation in our simulations of the subunit SDEs, although 
this approximation has been used in past simulation stud- 
ies [1, 5, 6]. We solve these systems of SDEs using the 
Euler-Maruyama method [41]. 



B. Simulation Results: Voltage Clamp 

In Figs. 2 and 3, we compare results from numerical 
simulations for a membrane patch size of 10 fim 2 against 
the analytical calculations presented in Section III. To 
simulate the voltage clamp condition, we fix V at a par- 
ticular value and keep it constant throughout the simu- 
lation. Figure 2 shows the mean and standard deviation 
of the proportion of open Na + and K + channels as a func- 
tion of the voltage clamp value. In most cases, the values 
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computed from numerical simulations (symbols) match 
the analytical results (lines). Of particular note is the 
fact that the computed values for the MC model (circle) 
and the channel model (x) are virtually indistinguishable. 

The only deviation between the numerical results and 
the analytical solutions occurs in the subunit models for 
the mean values of the Na + channels at high voltage val- 
ues. The cause of this discrepancy is that the analyt- 
ical treatment assumes that n is Gaussian-distributed 
whereas in the numerical methods the values of n are 
bounded between zero and one. For high voltage values 
of V, the proportion of open Na + channels is very small 
and the variance is non-zero so approximating the distri- 
bution of n as a Gaussian will allow n to take negative 
values. This cannot occur in the numerical simulations 
so the theoretical value for the mean fraction of open 
Na + channels will be less than the simulated value. As 
the number of channels increase, the variance of the frac- 
tion of open Na + channels decreases, which decreases the 
probability that a Gaussian-distributed n will take neg- 
ative values. The discrepancy between the analytically 
and numerically calculated values for the mean fraction of 
open Na + channels will decreases, therefore, as the num- 
ber of Na + channels increases. 



C. Simulation Results: Interspike Intervals 

Figure 4 shows mean and CV (left and right columns, 
respectively) of ISIs for three membrane areas that in- 
crease from top to bottom in each column. The input to 
the model is a constant DC input (I no i se — in Eqn. 26). 
The value of Idc is shown on the x-axis. In general, these 
simulations show that the rate and regularity of spiking 
activity produced by the MC (black line) and channel 
(gray line) models are in close quantitative agreement 
whereas the IdS (dashed line) and InS (dotted line) mod- 
els produce, on average, dramatically longer ISIs. The 
behavior of the models is most disparate for low stimu- 
lus levels and larger membrane areas. In these cases, the 
stimulus is not sufficiently strong to drive regular firing, 
so spike events are predominantly determined by stochas- 
tic fluctuations in the conductances due to channel noise. 
As the DC current is increased, the models respond to the 
external stimulus and there is a smaller effect of channel 
noise on spike timing. This leads all models to exhibit 
similar to mean ISI and CV values at high current levels. 

The mean ISIs for the subunit SDE models exceed 
those of the MC model for all stimulus conditions and 
membrane areas. This has been previously observed for 
the IdS model [6]. In the voltage clamp analysis we 
found that the variance in the proportion of open K + and 
Na + channels is much smaller for both of these models 
than for the MC model. This lack of conductance fluctu- 
ations leads to reduced spike rates at low stimulus levels, 
relative to the MC model. 

Overall, the results for the channel model show that 
it is possible to approximate the MC model with SDE 



and still obtain quantitatively accurate results, even for 
small numbers of channels. This is an important and a 
nontrivial result — the system size expansion is formally 
valid only in the limit of infinitely many channels, but 
we show here that it can be applied to a small number 
of channels, where membrane fluctuations have a major 
impact on spiking statistics. Moreover, the equilibrium 
noise approximation and treatment of boundary condi- 
tions do not appear to substantially degrade solution ac- 
curacy over a wide range of stimuli. 

Nevertheless, there are some discrepancies between the 
channel and MC models. At the smaller membrane areas 
(1 fim 2 and 10 /im 2 ), the mean ISIs for the channel model 
tend to be longer than the MC model ISIs. At the largest 
membrane area tested, for the case of weak or no input 
current, this trend is reversed and the channel model 
has shorter mean ISIs than the MC model. A possible 
source for these between the channel and MC models is 
our treatment of the boundary conditions in the channel 
model and the equilibrium noise approximation. Further 
investigation of this approximation is needed, but the 
similar ISI statistics obtained with the channel and MC 
models suggest that this approximation may be suitable 
in many cases. 

Figure 5 shows results obtained for two different levels 
of input noise added to the DC current amounts shown 
on the x-axis. We only present ISI statistics for the mem- 
brane size of 100 /im 2 because the smaller membrane ar- 
eas produce the same qualitative differences among the 
models. Overall, the effect of the stimulus noise is to re- 
duce the mean ISIs. Importantly, the ISI statistics of the 
responses of the MC and channel SDE models to these 
noisy stimuli remain quantitatively similar. In fact, the 
stimulus fluctuations elicit spikes even at low or no DC 
current levels, so the differences in the mean ISIs between 
the MC and channel SDE models become less apparent. 
Overall, this result indicates that the equilibrium noise 
approximation does not break down in the presence of a 
rapidly fluctuating external stimulus. Finally, the results 
for the subunit SDEs show, once again, that the stochas- 
tic dynamics and spiking activity of these approximate 
models do not accurately replicate the statistics of the 
MC model. 



V. DISCUSSION 

Beginning with the work of Fox and Lu [1, 29], the 
question of whether SDE models of channel noise can 
accurately approximate MC models has been explored. 
SDE models of membrane voltage fluctuations in HH 
models have several attractive features, including pos- 
sible improvements in the speed of numerical simulations 
and the opportunity to analyze these models using non- 
linear SDE theory [42, e.g.]. In recent years, however, 
the SDE approach has come under increasing scrutiny. 
Numerical simulations of the most commonly used SDE 
model, which we have called the identical subunit model, 
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FIG. 4. Means and coefficients of variation (CV) of first 2000 
interspike intervals as a function of the DC current level (ab- 
scissa) for a constant current input, (a) Results for a mem- 
brane area of 1 firn 2 . (b) Results for a membrane area of 
10 /im 2 . (c) Results for a membrane area of 100 fj.m 2 . 



have shown that this approach produces weaker conduc- 
tance and voltage fluctuations than the corresponding 
MC model [5, 31]. As a consequence, the firing rates of 
this SDE model are substantially lower [6] (and, equiva- 
lently, the mean ISIs are longer [3] ) , there is less variabil- 
ity in the occurrences and timing of spikes in response to 
a brief pulse of current [2] , and information is transmit- 
ted at a higher rate [6] . Furthermore, these discrepancies 
persist even as the number of channels increases [3, 5, 6]. 
In sum, there is an emerging consensus in the literature 
that the MC model cannot be approximated accurately 
using a subunit system of SDEs. 

We have demonstrated in this paper that an alternative 
SDE approach that is based on the multi-state structure 
of each ion channel can approximate the channel noise 
effects that are present in the MC model, even for rela- 
tively small numbers of channels, as long as the system- 
size expansion that is used to derive the SDE model is 
carried out properly. If one first defines the structure of 



FIG. 5. Means and coefficients of variation (CV) of first 2000 
interspike intervals as a function of the DC current level (ab- 
scissa) for a current input of the form Iuc + Inois^t where 
£t is a Gaussian white noise process with mean zero and unit 
variance. Membrane area is 100 fim 2 . (a) and (b) show 
results for I n0 i se = 1 and 2 fim 2 , respectively. 



a channel and then defines the dynamics of the propor- 
tions of channels in each configuration, one arrives at the 
channel-based SDE model (see discussion in Section II D 
and equations in Section IV A 2). If, instead, one approx- 
imates the proportion of subunits in the open or closed 
states with an SDE, then one obtains a subunit-based 
SDE model (see discussion in Section II D and equations 
in Section IV A 3) . Through our analysis of the stationary 
statistics of the proportion of open channels in voltage 
clamp, we have shown that the former approach, which 
we call channel-based, can provide a quantitatively accu- 
rate approximation to the MC model. We have also con- 
firmed that the latter, subunit-based approach, should 
not be considered an approximation of the MC model 
because it has fundamentally different stochastic prop- 
erties from the MC model. We conclude that the SDE 
approach is a valid approximate model for channel noise, 
but that is necessary to properly define the system of 
SDEs based on the structure of each channel. In particu- 
lar, one cannot include noise in the subunit equations in 
the manner suggested by Fox and Lu and expect results 
that are consistent with the MC ion channel model. 

To our knowledge, we are the first to present the nu- 
merical results of the channel SDE model [32] . We sim- 
ulated three membrane areas (1, 10, and 100 /im 2 as in 
[6]), where the number of Na + channels range from 60 to 
6000 and the number of K + channels range from 18 to 
1800. Our simulation results show that, in most cases, 
the ISI statistics for this model in response to constant 
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and noisy current inputs are in close quantitative agree- 
ment with the MC model (see Figs. 4 and 5). This finding 
is encouraging because the channel SDE model was de- 
rived by Fox and Lu using a system-size expansion [1, 29] 
that is formally valid only in the limit of a large number 
of channels. We did not find any evidence that the ap- 
proximation was invalid for the finite populations of ion 
channels at hand, and note that in many applications, 
the channel counts are in fact much higher. For instance, 
Rowat suggested that typical channel numbers in spike 
initiation zones may be on the order of 10 4 to 10 6 [17]. 
The channel counts used in the present study may be 
relevant to applications in which small patches of neu- 
ral membrane can drive spiking activity. For example, 
a node of Ranvier of the auditory nerve fiber can pro- 
duce a spike in response to cochlear implant stimulation. 
Typically, the nodes have surface areas of a few square 
micrometers [43] and are usually modeled with 1000 or 
fewer Na + channels [43-45]. 

As first pointed out by Fox and Lu [1] and as discussed 
in Section IV A 2, numerical simulations of the channel 
SDE model can be computationally expensive. One par- 
ticularly compute-intensive part of the algorithm is calcu- 
lating matrix square roots to determine stochastic terms 
in the SDE at every time step of the simulation. We have 
performed this operation using the optimized CBLAS li- 
brary, yet the channel SDE model still required approxi- 
mately 25 times as much computational time as the IdS 
model. Fortunately, as is the case with all three SDE 
approaches, the channel model has one considerable ad- 
vantage over MC — its computation time does not depend 
on the number channels. For example, in our implemen- 
tation of the Gillespie method, the computational time 
increased 12-fold as we increased the number of chan- 
nels from 600 Na+and 180 K+to 6000 Na+and 1800 K+. 
We found that even the channel model, slowest of the 
SDE approaches discussed in this paper, is faster than 
the MC model once the number of channels is greater 
than approximately 1200 Na + and 360 K + . Furthermore, 
the computational burden of solving the channel SDE 
model may be reduced by considering other methods for 
computing matrix square roots [46, e.g.]. Higher-order 
SDE solvers than the Euler-Maruyama method could also 
speed up SDE simulations [47, e.g.]. 

Even with increasingly efficient numerical methods for 
the channel SDE, a stochastic model that more directly 
resembles the classical HH equations, as opposed to the 
twelve-dimensional system of SDEs that defines the chan- 
nel model, is still desirable, as it would connect more 
directly with a wealth of studies of the original four- 
dimensional HH equations. In a simulation study of the 
Fox and Lu model, Bruce [5] sought to derive noise terms 
numerically for the subunit equations that would "cor- 
rect" the Fox and Lu model so that the fluctuations in 
the proportions of open channels would match the MC 
model. We have performed a similar analysis using the 
analytical voltage clamp results (see Appendix C for de- 
tails). We can redefine the magnitude of the noise in 



the subunit equations (i.e., a n in Eqn. 8) to produce a 
modified subunit SDE model with the same means and 
variances for the proportion of open channels in voltage 
clamp, to 0(1/ N). This model is constructed based on 
our analytical results for the stationary distributions of 
the proportions of open channels, so we call it a qua- 
sistationary model. The ISI statistics for this modified 
subunit SDE model are shown by the dashed line in Fig. 
6. This demonstrates that adjusting the noise terms in 
a subunit SDE model in order to fit the variances of the 
open channels in voltage clamp is not sufficient to pro- 
vide an improved fit to the spiking dynamics of the MC 
model. 

The reason for this can be seen in the multiple time 
scales of the autocorrelation functions for the MC model 
(Eqn. 22) and subunit models (Eqns. 24 and 25). We can 
alter the noise terms in the subunit model so that it pro- 
duces the correct stationary variances of the fractions of 
open channels, but we cannot modify the autocorrelation 
functions in a way that makes them consistent with the 
MC model. We therefore formulated a second quasista- 
tionary model by adding colored Gaussian noise to the 
conductances in the HH equations (details of this model 
arc in Appendix C). This second quasistationary model 
produces the stationary distribution of the channel SDE 
in voltage clamp, so the proportion of open channels in 
voltage clamp for this model has the same mean, vari- 
ance, and autocorrelation as the MC and channel SDE 
models. 

ISI statistics for this model are shown by the gray line 
in Fig. 6. We found that this model reproduced the 
statistics of the MC model much better than any of the 
subunit SDE models, so we conclude that temporal cor- 
relations in the channel noise play a critical role in in- 
fluencing spike timing. Moreover, since the structure of 
the ion channel determines the history-dependence of the 
channel noise, a valid channel noise model must prop- 
erly describe the dynamics of the entire channel and not 
only the kinetics of individual subunits. Using numeri- 
cal simulations, Bruce has pointed out that the subunit 
model does accurately approximate the MC model for 
the case of channels with a single subunit [5]. Our anal- 
ysis explains this observation because, for channels with 
one subunit, the channel-based and subunit-based SDE 
models are mathematically identical. 

We close by emphasizing that the models described 
in this paper do not represent complete descriptions of 
channel noise. As is always the case, when one attempts 
to formulate a mathematical descriptions of complex bi- 
ological processes, numerous assumptions and simplifica- 
tions are at play. As our understanding of the structure 
and dynamics of these membranes improves, it may be 
necessary to update and improve our mathematical mod- 
els of channel noise. Nonetheless, the central theme of 
this work will remain relevant: the statistics of channel 
noise are shaped by the activity of individual ion chan- 
nels, and therefore the approximation methods must also 
include information about the states of the channels in 
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FIG. 6. Means and coefficients of variation (CV) of first 2000 
interspike intervals for the MC model (black line) and two 
quasistationary models (see text and Appendix C for details) 
in response to constant current input. DC current level is 
given by the abscissa. Membrane area is 100 fxm 2 . 

order to correctly capture the state and history depen- 
dence of channel noise. 
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Appendix A: Matrices used in numerical simulations 
of the channel SDE model 

The state vectors are defined as x = [x\, x 2 , x 3 , Xi\ T 
and y = [y 1Q , y 2 o, J/30, yoi, Uu, 2/2 1, y3i] T ■ The number of 
K + and Na + channels arc denoted by Nk and A^a, re- 



spectively. The matrices Ak and in Eqn. 28 are: 



~(3a n +/3 n ) 
3a„ 





2/3„ 
-2(a n +/3„ 
2a„ 






3/3„ 

-(a„+3/3„) 4/3„ 

a„ -4/3„ J 





■ -(2*a m +P m +a h ) 


2,3 m 








Ph 










2a m 


-(a m +2/3 m +a h ) 


3l3 m 








0h 










&m 


~{3p m +a h ) 
























-(3a m +f) h ) 


m 










ce h 








3a m 


-(2a m +P m +P h ) 


2P m 


















2a m -{a r 


,+2/3 m +/3 h ) 


3/3 m 










a h 








Q!rn 


-(3/3 m +0 h )J 



The matrices Sk and SNa are the square root matrices 
of the following diffusion matrices: 



D* 



1 

Ak 



■4a„2 + (3a„+/3„)2i+2/3„S 2 -(3a„2i+2/3„2 2 ) 

-(3a„2i+2,3„S 2 ) 3a„2i+2(a„+/3„)2 2 +3/3„2 3 -(2a n x 2 +3/3„x 3 ) 

-(2a n x 2 +3/3 n x 3 ) 2a„S 2 + (a„+3/3„)s 3 +4^„24 -{a n x 3 +iP„Xi) 

— (q„x 3 +4,9„24) a n x 3 +i/3 n X4 . 



-2(a m y 10 +p m y 2a ) 



-("hJ/io+/3)ii/ii) 





-2(a m yio+0 m V2o) 
d 2 

-(a m y 20 +3P m y 30 ) 



-(a h y 20 +Phy 2 i) 




-(« 



V30 ) 





y 2 o+3/3. 
d 3 




(ahy 3 o+/3hy 3 i) 









-("hyio+^hj/n) o o 

-{a h y 20 +Phy 2 l) 

-(cthV '30+ Phil '31 ) 

-(3a m yoi+/3 m yn) 



-(3a m j/ 01 +/3 m |7i 1 ) d 5 -2(a m j7ii+/3 m 3/ 2 i) 

-2(a m yn+/3 m y 2 i) d 6 -(a m y 21 +3f3 m y 3 i) 

-(amS21+3/3 m i/3l) ^7 



where the elements on the diagonal are: 



di = 3a m j/ o + {2a r , 
d 2 = 2a m y w + (a m 
d3 = a m mo + (3/?m 
di = a h y o + (3a m 



+ m + a h )y w + 20 m y 2O + PhVu 

- 20 m + a h )y 20 + 3/3 m j/ 30 + h y 21 

- ah)V30 + PhV~3i 



d 5 = a h yw + 3a m j/oi + (2a m + m + 0h)yu + 20 m y 21 
d 6 = a h y 20 + 2a m y n + (a m + 2f3 m + f3 h )y 21 + 30 m y 31 
d 7 = a h y 30 + a m y 21 + (30 m + f3 h )y 31 

As discussed, we use the equilibrium mean values of x 
and y in the diffusion matrices. They are: 



x j 



4- 



<0 



K + p n ) 4 



and 



Vij = 



a 1 



3-i 



+ m ) 3 {a h + 0h)' 



Appendix B: Equilibrium Noise Approximation 

a. Voltage Clamp 

The equilibrium noise approximation in voltage clamp 
can be justified using a small noise expansion [34]. Let 
us define the small parameter e = 1/N and fix the mem- 
brane potential at a voltage clamp value so that a n and 
n can be treated as constants. Assume that n can be 
written as a series in the small noise parameter e: 

n(t) = n (t) + eni(i) + e 2 n 2 {t) ■ ■■ . 

If wc plug the small noise expansion into the subunit SDE 
for n in Eqn. 8 and collect terms of O(l), we find: 

driQ 
~dt 

and for terms of O(e): 
dn\ 



a n (l - no) - n n o , 



(Bl) 



dt 



= -a n m - n m + y/ a n (l - n ) + n n o £{t). (B2) 
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Equation Bl shows that hq satisfies a deterministic 
equation that does not depend on stochastic fluctuations 
in n. In the context of analyzing the stationary distribu- 
tion of n, therefore, we are justified in replacing no with 
its equilibrium value - . If we make this substitution 
for no in the equation for ri\ (Eqn. B2) and then form the 
sum no + eni, it is straightforward to arrive at the OU 
process in Eqn. 16. The equilibrium noise approximation 
is therefore the 0(1 /N) approximation of the long-time 
behavior of n(t) in voltage clamp. The same argument 
applies for the to and h subunits as well as for the mul- 
tivariate SDE that defines the channel SDE model. 



b. Time- dependent voltage 

Fox and Lu suggested applying this approximation in 
all cases, not just voltage clamp [1, 29] and we have 
used the approximation to simplify the numerical meth- 
ods for solving the channel SDE model. When V is not 
in voltage clamp, it evolves naturally and complicates 
the small noise expansion because it can introduce ad- 
ditional stochastic fluctuations into the gating variables 
and the voltage-dependent functions transition rate func- 
tions. Fox argued that the approximation would be ac- 
curate if the relaxation of V to its equilibrium value oc- 
curred on a much slower time scale than the relaxation 
of the gating variables (for the case of the subunit SDE 
model) [29] . Unfortunately, this separation of time scales 
does not appear to be a generic feature of HH models. 
Nevertheless, as shown in Figs. 4 and 5, the equilibrium 
noise approximation appears to be sufficiently accurate 
to reproduce spiking statistics to a high degree of accu- 
racy. 

Appendix C: Quasistationary Models 

In the Discussion (Section V), we introduced two mod- 
els that we discuss in greater detail here. We refer to 
both models as quasistationary approximations because 
they rely on results from our analysis of the station- 
ary statistics of open channels in Section III. In the 
first model, we were motivated by [5] to attempt to im- 
prove the accuracy of the subunit SDE model by modi- 
fying the noise terms in the gating equations. We have 
shown that the stationary variances for the proportion 
of open channels in the subunit SDE models does not 
match those of the MC model. To correct for this dis- 
crepancy, we can redefine <t^(V) in Eqn. 8 to guarantee 
that the stationary variance of n matches the stationary 
variance of the proportion of open K + channels under the 
MC models. The problem is simplified if we invoke the 
equilibrium noise approximation and use the facts that 
Var[n 4 ] = E[n 8 ] — E[n 4 ] 2 and that these higher moments 
of the stationary distribution of n are known since in 
voltage clamp n is an OU (Gaussian) process with mean 
A*sub = a a +a and variance a\ (V)t 71 (V)/2. The final 



step is to set Var[n 4 ] = Var[/ C h an ] and solve for a n (V). 
The exact solution would require inverting a nonlinear 
equation, but if we neglect terms that are of higher or- 
der than 0(1/N 2 ), we find a n (V) by finding the roots of 
168^ 4 o- 4 + l^Mn^n — — tAi)) I N , which is quadratic 

in er 2 . The same approach was also used to derive a new 
expression for a m so that this model also had the same 
stationary variance for the Na + channel. The formula for 
o h was left unchanged. 

As explained in the Discussion, this first approach 
did not provide a satisfactory approximation to the MC 
model, so we formulated a second quasistationary ap- 
proximation. We constructed this second model so that 
it and the MC model would have the same autocorrela- 
tion functions for the proportion of open channels and 
the same means and variances in voltage clamp. The 
mathematical structure of this model is somewhat un- 
usual in that the conductance is defined as the sum of 
a deterministic part (given by the HH equations for to, 
h, and n) and a colored Gaussian processes that is de- 
fined by the autocovariance function for the proportion of 
open channels in the MC model. As usual, we illustrate 
our approach with the K + channel. The conductance is 
defined to be: 



9k = .9K(n 4 + o- chan 7] t ) (CI) 



where n is the classical (deterministic) gating variable 
satisfying an ordinary differential equation of the form 
of Eqn. 3 and r\ t is a stochastic process. To define rjt, 
first note that the channel SDE model provided a quan- 
titatively accurate approximation of the MC model, so 
it is reasonable to describe the stationary conductance 
as a Gaussian process. Second, recall from Eqn. 22 that 
the autocorrelation function for the proportion of open 
K + channels has four distinct time scales (the first four 
multiples of l/r„). Taken together, these facts lead us 
to model the K + conductance as a non-Markovian Gaus- 
sian process. The representation theory of Gaussian pro- 
cesses furnishes a systematic method for constructing the 
stochastic process % based on the autocorrelation func- 
tion for the K + channel [48]. In particular, this theory 
guarantees that r\ t can be written in terms of a stochas- 
tic (Wiener) integral of the form: 



f 4 i(t-u) 

J o i=l 



The coefficients dj are voltage-dependent and are com- 
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puted by solving the system of nonlinear equations: 
r n (V)a 1 (V)J2 a ^l=4a n (Vfr 

i=l 

T n (V)a 2 (V) V = 6a n (V) 2 f3 n (V)T 

1— l 

r n (V)a 3 (V)J2 = 4a„(U)/3„(U) 2 r 

2— 1 

r n (T/)a 4 (y)^^=/?„(t/) 3 r 

i=l 

where T = {an{ y)+[ij!v)y- an {vy ■ In Practice. the sys- 
tem of nonlinear equations defining dj(U) must be solved 
numerically. We used the Minimize command in Maple 
(Waterloo Maple Inc., Version 13) to generate a data 
table of values of (ii(V) that solved these equations for 
voltage values ranging from —20 to 120 mV in increments 
of 0.01 mV. The procedure for constructing the non- 
Markovian Gaussian process for the Na + conductance is 



similar but slightly more complicated because there are 
seven time scales in the autocorrelation function. We 
omit these details here and direct the interested reader 
to computer code available at ModelDB [35]. 

To numerically integrate this quasistationary model, 
we used a forward Euler method to update the value of 
V and n in each time step, where the stochastic integral 
for rjt is integrated as follows: 

1. Compute the voltage dependent terms r„, <r n , a n , 
/?„, and oij using the value of V from the previous 
time step. 

2. Update the terms associated with each time scale: 
Ai(t + At) = Ai(t)e~~ + aiG n r\J~K~t where i = 
1, 2, 3, 4, At is the time step, and r is a mean zero, 
unit variance Gaussian random generated on each 
time step. 

3. Update the stochastic process: = 
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